classdef CELM_model < handle
    % MELM_MODEL 改进的极限学习机（Modified Extreme Learning Machine）模型
    %
    % 此类实现了MELM模型的训练和测试功能，专为分类任务设计。
    % 特点：
    %   - 使用改进的权重初始化策略（基于类间样本关系）
    %   - 支持多种激活函数（sigmoid, sine, hardlim, tribas, radbas）
    %   - 自动处理多分类标签编码
    %   - 计算训练/测试时间和准确率
    
    properties
        Task_Type                   % 任务类型标识 (0:回归, 1:分类)
        NumberofHiddenNeurons       % 隐含层神经元数量
        ActivationFunction          % 激活函数类型（字符串）
        C                           % 正则化参数（控制模型复杂度）
        Solver_flag                 % 求解器标志位（0:伪逆求解, 1:QR分解求解）
    end
    
    properties (Access = private)
        InputWeight                 % 输入层到隐含层的权重矩阵 [hidden_neurons x features]
        BiasofHiddenNeurons        % 隐含层神经元的偏置向量 [hidden_neurons x 1]
        OutputWeight               % 隐含层到输出层的权重矩阵 [output_neurons x hidden_neurons]
        label                      % 训练集的唯一标签（用于测试时编码）
        train_time_                % 训练耗时（秒）私有存储
        test_time_                 % 测试耗时（秒）私有存储
        train_accuracy_            % 训练精度私有存储（0-1范围）
        test_accuracy_             % 测试精度私有存储（0-1范围）
    end
    
    properties (Dependent)
        TrainingTime                % 获取训练耗时（秒）的依赖属性
        TestingTime                 % 获取测试耗时（秒）的依赖属性
        TrainingAccuracy            % 获取训练精度的依赖属性
        TestingAccuracy             % 获取测试精度的依赖属性
    end
    
    methods
        function obj = CELM_model(Task_Type, NumberofHiddenNeurons, ActivationFunction, C, Solver_flag)
            % MELM_MODEL 构造函数初始化MELM模型参数
            %
            % 输入参数：
            %   Task_Type             - 任务类型 (0:回归, 1:分类)
            %   NumberofHiddenNeurons - 隐含层神经元数量（正整数）
            %   ActivationFunction    - 激活函数名称（'sigmoid','sine','hardlim','tribas','radbas'）
            %   C                     - 正则化参数（标量，值越大正则化越弱）
            %   Solver_flag           - 求解器选择（0:伪逆, 1:QR分解）
            
            % 初始化模型参数
            obj.Task_Type = Task_Type;
            obj.NumberofHiddenNeurons = NumberofHiddenNeurons;
            obj.ActivationFunction = ActivationFunction;
            obj.C = C;
            obj.Solver_flag = Solver_flag;
        end
        
        function train(obj, train_data)
            % TRAIN 训练MELM模型（使用改进的权重初始化策略）
            %
            % 输入参数：
            %   train_data - 训练数据矩阵 [样本数 x (标签+特征)]，第一列为标签
            
            % 常量定义（回归任务标识）
            REGRESSION = 0;

            % 分离标签和特征
            T = train_data(:,1)';              % 标签行向量 [1 x samples]
            P = train_data(:,2:end)';           % 特征矩阵 [features x samples]
            
            % 分类任务：预处理标签（one-hot编码并转换为-1/1格式）
            if obj.Task_Type ~= REGRESSION
                % 调用私有方法进行标签编码
                [T, obj.label] = obj.preprocess_labels(T);
            end

            % 改进的权重初始化策略（获取输入权重和隐含层偏置）
            [obj.InputWeight, obj.BiasofHiddenNeurons] = obj.init_weights(P, T);
            
            % 计算隐含层输出矩阵 H
            start_time = cputime;  % 开始计时
            
            % 计算加权输入：权重矩阵 * 特征矩阵
            tempH = obj.InputWeight * P;
            % 扩展偏置向量为矩阵（与样本数匹配）
            BiasMatrix = repmat(obj.BiasofHiddenNeurons, 1, size(P,2));
            % 添加偏置
            tempH = tempH + BiasMatrix;
            % 应用激活函数得到隐含层输出 H [hidden_neurons x samples]
            H = obj.activate(tempH, obj.ActivationFunction);
            clear P tempH;  % 清理中间变量释放内存
            
            % ========= 计算输出权重 =========
            if obj.C == 10^100  % 无正则化情况
                if obj.Solver_flag == 0
                    % 伪逆求解：OutputWeight = pinv(H') * T'
                    obj.OutputWeight = pinv(H') * T';
                else
                    % QR分解求解（调用外部函数）
                    obj.OutputWeight = QRIL_inverse.QR_schmidt(H', T');
                end
            else  % 带正则化求解
                if obj.Solver_flag == 0
                    % 正则化伪逆求解：(I/C + H*H')^{-1} * (H * T')
                    obj.OutputWeight = (eye(size(H,1)) / obj.C + H * H') \ (H * T');
                else
                    % QR分解求解正则化系统
                    obj.OutputWeight = QRIL_inverse.QR_schmidt(...
                        (eye(size(H,1)) / obj.C + H * H'), eye(size(H,1))) * H * T';
                end   
            end
            
            % 记录训练结果
            obj.train_time_ = cputime - start_time;  % 计算训练耗时
            
            % 计算训练集预测输出：Y = (H' * OutputWeight)'
            Y = (H' * obj.OutputWeight)';
            % 计算训练准确率
            obj.train_accuracy_ = obj.calculate_accuracy(T, Y);
        end
        
        function test(obj, test_data)
            % TEST 测试MELM模型性能
            %
            % 输入参数：
            %   test_data - 测试数据矩阵 [样本数 x (标签+特征)]，第一列为标签
            
            % 分离标签和特征
            TV_T = test_data(:,1)';     % 测试标签
            TV_P = test_data(:,2:end)';  % 测试特征
            
            % 获取测试样本数量
            NumberofTestingData = size(TV_P, 2);
            
            % 编码测试标签（使用训练时保存的标签映射）
            TV_T = obj.encode_labels(TV_T);
            
            % 计算测试集隐含层输出
            start_time = cputime;  % 开始计时
            
            % 计算加权输入
            tempH_test = obj.InputWeight * TV_P;
            % 扩展偏置向量
            BiasMatrix = repmat(obj.BiasofHiddenNeurons, 1, NumberofTestingData);
            % 添加偏置
            tempH_test = tempH_test + BiasMatrix;
            % 应用激活函数
            H_test = obj.activate(tempH_test, obj.ActivationFunction);
            clear TV_P tempH_test;  % 清理中间变量
            
            % 计算预测输出：TY = (H_test' * OutputWeight)'
            TY = (H_test' * obj.OutputWeight)';
            
            % 记录测试结果
            obj.test_time_ = cputime - start_time;  % 计算测试耗时
            % 计算测试准确率
            obj.test_accuracy_ = obj.calculate_accuracy(TV_T, TY);
        end
        
        %% 依赖属性的访问方法（getter）
        function value = get.TrainingTime(obj)
            % 获取训练耗时
            value = obj.train_time_;
        end
        
        function value = get.TestingTime(obj)
            % 获取测试耗时
            value = obj.test_time_;
        end
        
        function value = get.TrainingAccuracy(obj)
            % 获取训练精度
            value = obj.train_accuracy_;
        end
        
        function value = get.TestingAccuracy(obj)
            % 获取测试精度
            value = obj.test_accuracy_;
        end
    end
    
    methods (Access = private)
        function [InputWeight, BiasofHiddenNeurons] = init_weights(obj, P, T)
            % INIT_WEIGHTS 改进的权重初始化策略
            %
            % 输入参数：
            %   P - 特征矩阵 [features x samples]
            %   T - 编码标签 [classes x samples]（分类任务）
            %
            % 输出：
            %   InputWeight - 初始化后的输入权重矩阵
            %   BiasofHiddenNeurons - 初始化后的隐含层偏置向量
            
            % 初始化输出变量
            InputWeight = zeros(obj.NumberofHiddenNeurons, size(P,1));
            BiasofHiddenNeurons = rand(obj.NumberofHiddenNeurons, 1);
            
            % 获取样本数量和类别数量
            [~, N] = size(P);  % N = 样本数
            [number_class, ~] = size(T);  % 类别数（仅分类任务有效）
            
            % 定义常量
            REGRESSION = 0;  % 回归任务标识
            Flag = 1;        % 权重过滤控制变量
            
            % ===== 分类任务权重初始化 =====
            if obj.Task_Type ~= REGRESSION 
                i = 1;
                % 循环初始化每个隐含层神经元
                while i <= obj.NumberofHiddenNeurons
                    
                    % 异类样本差 
                    % 获取异类样本索引
                    same_class = false;
                    [idx1, idx2] = obj.getValidSamples(T, same_class);
                    
                    % 计算权重向量（异类样本特征差）
                    w = P(:, idx1) - P(:, idx2);
                    norm_w = dot(w, w);  % 计算范数平方
                    
                    % 权重过滤：跳过过大的权重
                    if norm_w < 1/(Flag * number_class)
                        Flag = Flag + 1;
                        continue;
                    end
                    
                    % 归一化并存储权重
                    InputWeight(i, :) = w / (norm_w/2);
                    
                    % 计算特殊偏置：基于样本范数差
                    BiasofHiddenNeurons(i,1) = ...
                        (dot(P(:,idx2), P(:,idx2)) - dot(P(:,idx1), P(:,idx1))) / norm_w;
                  
                    i = i + 1;  % 成功初始化一个神经元，索引+1
                end

            % ===== 回归任务权重初始化 =====
            else
                i = 1;
                while i <= obj.NumberofHiddenNeurons
                    % 随机选择两个不同样本
                    idx1 = randi(N);
                    idx2 = randi(N);
                    if idx1 == idx2
                        continue;  % 跳过相同样本
                    end
                    
                    % 计算样本特征差
                    w = P(:, idx1) - P(:, idx2);
                    norm_w = dot(w, w);  % 范数平方
                    
                    % 权重过滤
                    if norm_w < 1/(Flag + 10)
                        Flag = Flag + 1;
                        continue;
                    end
                    
                    % 归一化并存储权重
                    InputWeight(i, :) = w / (norm_w/2);
                    
                    % 计算偏置：基于样本范数差
                    BiasofHiddenNeurons(i,1) = ...
                        (dot(P(:,idx2), P(:,idx2)) - dot(P(:,idx1), P(:,idx1))) / norm_w;
                    
                    i = i + 1;  % 成功初始化一个神经元
                end
            end
        end

        function [idx1, idx2] = getValidSamples(obj, T, sameClass)
            % GETVALIDSAMPLES 获取满足类别条件的样本对
            %
            % 输入：
            %   T        - 编码标签矩阵 [classes x samples]
            %   sameClass - 布尔值：true=同类样本，false=异类样本
            %
            % 输出：
            %   idx1, idx2 - 满足条件的样本索引

            [~, N] = size(T);  % 获取样本数量
            attempts = 0;
            max_attempts = 1000;  % 防止无限循环
            
            % 循环尝试寻找有效样本对
            while attempts < max_attempts
                idx1 = randi(N);
                idx2 = randi(N);
                
                % 跳过相同样本
                if idx1 == idx2
                    continue;
                end
                
                % 获取样本的真实类别
                [~, cls1] = max(T(:, idx1));  % 最大概率对应的类别
                [~, cls2] = max(T(:, idx2));
                
                % 检查是否满足同类/异类条件
                if (sameClass && (cls1 == cls2)) || (~sameClass && (cls1 ~= cls2))
                    return;  % 找到有效样本对，直接返回
                end
                
                attempts = attempts + 1;  % 增加尝试计数
            end
            
            % 达到最大尝试次数：随机返回不同样本
            idx1 = randi(N);
            idx2 = randi(N);
            while idx1 == idx2
                idx2 = randi(N);  % 确保两个索引不同
            end
        end
        
        function H = activate(~, tempH, ActivationFunction)
            % ACTIVATE 应用激活函数
            %
            % 输入：
            %   tempH - 加权输入矩阵
            %   ActivationFunction - 激活函数名称字符串
            %
            % 输出：
            %   H - 激活后的输出矩阵
            
            switch lower(ActivationFunction)
                case {'sig','sigmoid'}
                    H = 1 ./ (1 + exp(-tempH));  % Sigmoid激活
                case {'sin','sine'}
                    H = sin(tempH);  % 正弦激活
                case 'hardlim'
                    H = double(tempH >= 0);  % 硬极限激活
                case 'tribas'
                    H = max(1 - abs(tempH), 0);  % 三角基函数
                case 'radbas'
                    H = exp(-(tempH.^2));  % 径向基函数
                otherwise
                    error('不支持的激活函数: %s', ActivationFunction);
            end
        end
        
        function [encoded_labels, unique_labels] = preprocess_labels(~, raw_labels)
            % PREPROCESS_LABELS 预处理分类标签（训练阶段）
            %
            % 输入：
            %   raw_labels - 原始标签向量 [1 x samples]
            %
            % 输出：
            %   encoded_labels - 编码后的标签矩阵 [classes x samples]（-1/1格式）
            %   unique_labels - 唯一的标签类别
            
            % 获取唯一标签和类别数量
            unique_labels = unique(raw_labels);
            num_classes = length(unique_labels);
            
            % 初始化编码矩阵
            encoded_labels = zeros(num_classes, length(raw_labels));
            
            % 构建one-hot编码
            for i = 1:length(raw_labels)
                class_idx = (unique_labels == raw_labels(i));  % 找到标签位置
                encoded_labels(class_idx, i) = 1;  % 设置对应位置为1
            end
            
            % 将[0,1]编码转换为[-1,1]格式
            encoded_labels = encoded_labels * 2 - 1;
        end
        
        function encoded_labels = encode_labels(obj, raw_labels)
            % ENCODE_LABELS 编码测试标签（测试阶段）
            %
            % 输入：
            %   raw_labels - 原始测试标签向量 [1 x samples]
            %
            % 输出：
            %   encoded_labels - 编码后的标签矩阵 [classes x samples]（-1/1格式）
            
            % 初始化编码矩阵（使用训练时保存的类别信息）
            encoded_labels = zeros(length(obj.label), length(raw_labels));
            
            for i = 1:length(raw_labels)
                % 查找标签在训练集中的位置
                class_idx = (obj.label == raw_labels(i));
                
                % 检查是否为未知标签
                if ~any(class_idx)
                    error('测试标签 %d 未在训练集中出现', raw_labels(i));
                end
                
                % 设置对应位置为1
                encoded_labels(class_idx, i) = 1;
            end
            
            % 转换为[-1,1]格式
            encoded_labels = encoded_labels * 2 - 1;
        end
        
        function accuracy = calculate_accuracy(~, T, Y)
            % CALCULATE_ACCURACY 计算分类准确率
            %
            % 输入：
            %   T - 真实标签矩阵 [classes x samples]
            %   Y - 预测输出矩阵 [classes x samples]
            %
            % 输出：
            %   accuracy - 分类准确率（0-1）
            
            miss = 0;  % 错误计数
            
            % 遍历所有样本
            for i = 1:size(T, 2)
                % 获取真实类别索引
                [~, expected] = max(T(:,i));
                % 获取预测类别索引
                [~, actual] = max(Y(:,i));
                
                % 统计错误样本
                if actual ~= expected
                    miss = miss + 1;
                end
            end
            
            % 计算准确率
            accuracy = 1 - miss / size(T, 2);
        end
    end
end